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The  finite  element  method  is  extended  to  the  refined  elastic  analysis 
of  multilayer  beams  and  shells  with  no  restriction  placed  upon  the 
ratios  of  the  layer  thicknesses  and  properties.  The  method  is  appli¬ 
cable  to  structures  wherein  shearing  deformations  are  significant, 
including  sandwich-type  structures.  Element  stiffnesses  developed 
are  based  on  polynomial  displacement  models  and  are  applicable  to 
the  linear  elastic  analysis  of  beams  and  thin,  axisymmetric  shells 
of  arbitrary  meridian.  Although  attention  is  restricted  to  three¬ 
layered  construction  with  similar  facings,  the  theory  may  be  gener¬ 
alized  for  any  type  of  flexural  element  and  any  arrangement  of 
laminations.  Computer  programs  have  been  written  for  both  static 
analysis  and  free  vibration  analysis.  Inclusion  of  rotary  as  well  as 
translational  inertia  allows  determination  of  natural  thickness -shear 
frequencies  and  mode  shapes  in  addition  to  flexural  vibration  charac¬ 
teristics.  Examples  are  presented  to  illustrate  the  effectiveness  of 
the  method. 
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SECTION  I 
INTRODUCTION 


Multilayer  construction  has  become  more  and  more  important  in  structural  engineering 
as  one  means  of  achieving  a  beneficial  combination  of  the  properties  of  two  or  more  materials. 
The  best  known  examples  of  this  type  are  the  widespread  “sandwich”  structures  used  in  the 
aerospace  industry.  These  combine  thin,  high-strength  facing  layers  with  a  thicker,  light¬ 
weight  core. 

The  theory  of  stress  analysis  of  multilayer  structures  is  well  established.  In  general, 
two  classifications  of  such  structures  can  be  identified:  (1)  “laminates”  in  which  layers 
of  materials  with  similar  properties  are  bonded  together  and  for  which  the  Kirchhoff-Love 
hypothesis  is  applied  and  (2)  “sandwiches”  in  which  some  layers  may  be  significantly  weaker 
than  others  and  for  which  transverse  shear  deformation  is  taken  into  account.  Theory  for  the 
laminates  (Reference  1)  has  been  successfully  applied  to  analysis  of  general  plates  and  shells 
using,  for  example,  the  approximate  methods  of  finite  differences  (Reference  2)  and  finite 
elements  (Reference  3).  However,  despite  the  availability  of  sandwich  theories  of  various 
degrees  of  refinement  in  the  literature,  there  have  been  relatively  few  solutions  published 
that  include  the  effects  of  transverse  shear  (Reference  23).  Moreover,  these  solutions  have 
been  restricted  to  the  simpler  geometries  such  as  rectangular  and  circular  plates  and 
cylindrical  and  spherical  shells. 

The  purpose  of  this  present  paper  is  to  extend  the  finite  element  method  to  the  analysis 
of  sandwich  beams  and  shells.  Although  specific  solutions  are  to  be  presented  here  only  for 
axisymmetric  shells,  the  same  general  approach,  using  existing  standard  finite  element 
techniques,  will  permit  the  analysis  of  arbitrary  configurations  and  boundary  conditions. 

Extensive  reviews  and  bibliographies  of  the  theory  of  sandwich  structures  are  presented 
in  References  4  through  6.  In  the  present  work,  a  theory  analagous  to  Yu’s  sandwich  theory 
is  adopted  (Reference  7).  This  theory  places  no  restriction  on  the  ratios  of  layer  thicknesses 
and  material  properties,  and  it  includes  the  bending  and  stretching  effects  and  the  transverse 
shear  flexibility  of  all  layers.  Yu  has  applied  this  formulation  to  vibrational  problems  of 
sandwich  structures,  including  both  shear  and  rotatory  inertia  (References  8  and  9). 

The  finite  element  method  was  developed  concurrently  with  the  increasing  use  of  high¬ 
speed  electronic  digital  computation  and  its  concommitant  emphasis  on  discretized  techniques 
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in  structural  analysis.  In  general,  the  method  has  proved  to  be  a  successful  tool  for  the 
systematic  analysis  of  complex  structures  and  the  approximate  solution  of  difficult  problems 
in  continuum  mechanics.  Among  the  large  number  of  publications  on  the  finite  element  method 
to  appear  in  the  last  decade,  References  10  and  11  include  descriptions  of  the  basic  techniques 
and  comprehensive  bibliographies. 

Because  of  the  importance  of  the  problem,  considerable  attention  has  been  devoted  to  the 
finite  element  analysis  of  axisymmetric  shells.  This  class  of  problems  has  the  inherent 
simplification  of  a  one -dimensional  mesh.  For  the  analysis  the  conical  frustra  as  elements 
have  been  widely  used  (References  12  through  14).  Recently,  the  more  accurate  axisymmetric 
types  of  doubly  curved  elements  have  been  introduced  (References  15  through  17).  In  this 
paper,  the  doubly  curved  element  developed  by  Khojasteh-Bakht  (Reference  17)  is  utilized. 
This  element  duplicates  the  position,  slope,  and  curvature  of  the  original  structure  at  each 
nodal  circle.  In  addition,  this  element  appears  to  be  satisfactory  from  the  point  of  view  of 
the  Irons-Draper  conditions  of  completeness  and  compatibility  (Reference  18).  In  this  devel¬ 
opment,  polynomial  displacements  are  assumed  and  the  direct  stiffness  method  is  employed. 

The  following  assumptions  apply  throughout  this  paper: 

1.  Displacements  and  strains  are  sufficiently  small  so  that  the  linear  theory  of  elasticity 
applies. 

2.  Perfect  bonding  occurs  between  adjacent  layers  of  the  structure. 

3.  The  transverse  displacement  of  all  layers  is  the  same  at  a  given  location  of  the 
middle  surface  of  the  structure.  In  other  words,  there  is  no  pinching  deformation. 

4.  Shells  are  thin  in  the  sense  that  products  of  thickness  with  curvature  are  much 
smaller  than  unity  (£  /R  «1). 

5.  Material  lines  in  each  layer  originally  straight  and  normal  to  the  neutral  surface 
remain  straight  after  deformation,  but  no  longer  remain  normal.  The  difference  in  shear 
strain  in  the  several  layers  manifests  itself  in  warping  of  the  cross  section  at  the  interfaces. 

6.  The  materials  of  each  layer  are  linearly  elastic  and  the  procedure  can  be  applied  to 
orthotropic  materials.  However,  for  simplicity,  the  discussion  in  this  paper  is  confined  to  an 
isotropic  case. 
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7.  The  bending  and  stretching  of  all  layers  are  taken  into  account.  However,  this 
assumption  can  be  relaxed  for  a  particular  layer  by  assigning  a  zero  Young’s  modulus. 

8.  All  layers  are  flexible  in  shear  (see  item  5  above)  but  this  assumption  can  be  relaxed 
for  a  particular  layer  by  assigning  an  infinite  shear  modulus. 


SECTION  n 

SANDWICH  BEAM  ELEMENT 

A  sandwich  beam  element  of  unit  width  which  is  symmetric  about  the  middle  surface  of 
the  core  is  shown  in  Figure  1.  A  normalized  coordinate  is  defined 

£  =  (x-x.)/(x.-x.)  =  (x  —  x,  )/jt 

where  the  subscripts  identify  the  ith  and  jth  nodes  of  the  element  of  length  £ .  The  slope  of 
beam  can  be  divided  into  contributions  due  to  pure  bending  deformation  and  pure  shearing 
deformation  as  follows: 

”df"  =  X  =  X5  X  s  (I) 

Under  the  above  assumptions  the  shear  strain  of  each  layer  is  independent  of  the  thickness 
coordinate  z.  The  shear  strains  of  the  faces  and  core  are  indicated  by  y^  and  yQ,  respectively. 
Figure  2  shows  a  differential  element  deforming  under  pure  shearing  of  the  core  and  the 
facings.  With  this  type  of  deformation  there  is  no  net  extension  of  the  layers,  so  the  axial 
displacement  of  each  middle  surface  is  zero.  The  axial  displacement  of  the  interface  must  be 
the  same  when  computed  with  reference  to  the  middle  surface  of  either  the  face  or  core.  For 
the  case  of  constant  shear  carried  entirely  by  the  core,  this  condition  gives 

u(z  =  h  _/2)  -(  y  ~  X  )hc/2  :  y  hf/2 
c  'c  Asc  c  A  s  c  * 


216 


AFFDL-TR-68-150 


and  for  the  case  of  face  shearing  it  gives 

u  (  z  =  -hc/2)  =  <  yt  ~Xsf)hf  /2  =  Xsf  hc/Z 

Adding  the  two  equations  and  using  d  =  h  +  h„,  one  obtains 

C  I 

xs  1  Xsc  +  xsf  =  xchc/d4/f  hf /d  (2) 

x  s£  and  xsc  are  defined  in  Figure  2. 


As  a  result  of  the  kinematic  assumptions,  the  axial  and  transverse  displacements  of  the 
beam  may  be  written 

uc  a  “  zc  (  X  ~  Xc  ) 


t ,b  ,  .  .  d  _  hc  _  hf 

'f  :  “zf  (x~xf  iJix+yyc  +  TXf 

wc  =  wf  =  wfb  =  w 


(3) 


where  the  subscripts  c  and  f  indicate  core  and  facings  and  the  superscripts  t  and  b  indicate 
the  top  and  bottom  facings.  zq  and  zf  are  thickness  coordinates  with  origins  at  the  middle 
surface  of  the  respective  layers.  The  strain-displacement  equations  for  the  beam  are 

duc 


"xc 


.t,b 
;  xf 


Xz 


xc 


t.b 

7 

7  z  x  f 


dx 

€. 

.  ».b 

auf 

dx 

duc 

d  w 

dz 

dx 

du* 

=  — -  +■ 

dw 

o 

xc 


+  zc  *xc 


■t,b 


dz 


xf 


=  X. 


dx  1  yf 


4-  7  1  ^ 

^  zfKxi 


(4) 


Substituting  Equations  3  into  Equations  4,  the  strain  components  that  result  are  given  by 


€  *  0 
xc 


,o  t,b 

.  +  d  dX 

5  JV 

xf 

'  "  2  dx 

2 

K 

*1 

1 

II 

irce 

xc 

dx 

dx 

t.b 

dX 

dXf 

K  ' 
xf 

~d7  + 

d  x 

i  +  2  d  x 


(5) 
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Figure  1.  Typical  Beam  Element 


Figure  2.  Slope  Due  to  Shearing 
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Finally,  the  stress-strain  equations  in  matrix  form  are 

r  cr 

xc 

*  Ec 

0 

0 

0 

0 

o' 

C 

*xc 

r 

zxc 

0 

0 

0 

0 

0 

^zxc 

cr  :  - 

7  : 

0 

0 

Ef 

0 

0 

0 

4 

t 

€xf  l 

Tzxf 

0 

0 

0 

\fGf 

0 

0 

4f 

crb 
x  f 

0 

0 

0 

0 

Ef 

0 

«b. 

xf 

rb 

.  zxf  J 

0 

0 

0 

0 

0 

XfGf 

y  b 
/zxf 

Ce 


(6) 


Here  \c  and  \f  are  shear-stress  correction  factors  analagous  to  that  used  in  the  Timoshenko 
beam  theory.  Yu  (Reference  7)  has  suggested  values  very  close  to  unity  for  sandwich  con¬ 
struction  with  thin,  heavy  facings  and  a  light,  weak  core.  Since  these  are  the  properties 
characteristic  of  the  vast  majority  of  sandwich  structures,  these  correction  factors  will 
hereafter  be  taken  as  unity.  However,  it  should  be  noted  that  in  a  generalization  of  the  present 
approach,  further  investigations  are  necessary  to  evaluate  these  parameters. 


The  beam  element  stiffness  matrix  is  derived  using  a  cubic  transverse  displacement  field 
and  a  linear  variation  of  shear  rotation.  Moreover,  interpolation  functions  are  used  in  order 
to  express  the  displacement  models  directly  in  terms  of  the  nodal  displacements.  Hence 


u(  £  )  =  <J>  (£)  q  ,  0<  £  <  l 


Here  the  vectors  are  chosen  as 


and  thus  the  matrix  <J>  is  given  by 


<  *  X  Yc  Yf  > 


<  u  (0) 


u  (I  )  > 


{7) 


( l-3f2+2£3) 
&£  (£-\)/£ 
o 
o 


-4fd-2£+£2) 

<l-4£+3£2  ) 
0 
0 


0  o  j  C  <3-20  i£2<£-l  ) 

0  0  |  G&\-£)/£  £<3£~2) 


<1-0  0  I  0  0 

I 

0  ( I  -  c  > 1  0  0 


0 

0 

£ 

0 


0 

0 

0 


£ 


By  applying  Equations  4  and  5  to  Equation  7,  the  strains  may  be  expressed  in  terms  of  the 
nodal  displacements 


«(  £  )  =  B  (  £  )  q 


<  8) 
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where  the  strain-component  vector  is  defined 


<  K 


xc  Xxzc 


€xf 


Kxf 


Xxzf 


rxf  *xf 


Xxbzf 


Then  consistent  with  the  definitions  in  Equation  6,  it  is  possible  to  write 


zc 

0  0  0 

0  0 

0 

o’" 

0 

1  0  0 

0  0 

0 

0 

e  =  Z  e  = 

0 

0  1  Zf 

0  0 

0 

0 

€ 

0 

0  0  0 

1  0 

0 

0 

0 

0  0  0 

0  1 

zf 

0 

0 

0  0  0 

0  0 

0 

The  strain  energy  of  the  element  is  now  given  by 

,  h/2 

»-U*T 

cr  dV  = 

^4T/  /  B(C)TZ( 

z )  CZ(z)B(£)dz 

d£  q 

0  _h/2 

Therefore,  the  element  stiffness  with  reference  to  the  displacements  4  is 


1  z 

k q  =  J  B(£)Tf  Z(  z  )T  CZ  (z  )dz  B(£  )  (10) 

8x8  0  _t\/z 

The  inner  integral  is  readily  recognized  as  a  diagonal  matrix  involving  the  material  proper¬ 
ties  and  the  areas  and  moments  of  inertias  of  the  layers.  The  8x8  stiffness  matrix  thus  has 
five  distinct  components  involving  core  bending  and  shear  and  facing  bending,  shear  and 
extension.  These  are  identified  by  the  E  Ic»  AcGc>  E^,  AfGf  and  AfEf  terms,  respectively. 


In  assembling  finite  elements  into  a  representation  of  an  overall  structure,  compatibility 
usually  must  be  maintained  for  all  the  displacement  degrees  of  freedom  occurring  at  the 
interelement  nodes.  However,  when  the  transverse  shear  behavior  is  included,  some  continuity 
conditions  must  be  removed  in  order  to  permit  the  “kinking”  associated  with  discontinuities 
of  the  shear  stress  resultant.  These  discontinuities  occur  at  transverse  concentrated  loads. 
Thus,  the  necessary  and  sufficient  requirement  for  compatibility  of  the  assemblage  is 
interelement  continuity  on  only  the  following  nodal  displacements: 

1.  the  transverse  displacement,  w 

2.  the  rotations  of  the  normals  to  the  reference  surface  (i.  e.,  the  rotations  associated 
with  pure  bending) ,  X  b 
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3.  the  shear  warping  angles  at  the  interfaces,  Y  -  Y  - 

C  I 

Hence,  the  stiffness  k^  must  be  transformed  so  that  it  is  expressed  in  terms  of  the  following 
nodal  displacements 

rT  =  '  "i  xw  y,  «,  xbj  y, |yf ,  yfj  >  mi 

The  transformation  is  quite  simple  and  can  be  constructed  from  the  definitions 

x  1  ye  "  * 

xb  =  x-hcrc  /<■-*•,  rf  /«• 

The  latter  is  obtained  from  Equation  2.  A  static  condensation  is  performed  on  the  resulting 
stiffness  to  obtain  the  6x6  portion  that  is  employed  in  the  direct  stiffness  method. 

It  should  be  added  that  the  nodal  displacement  parameters  given  in  Equation  11  are  also 
more  advantageous  than  those  of  Equation  7  for  expressing  the  boundary  conditions  necessary 
at  a  fixed  support.  At  such  a  location,  the  transverse  displacement  w,  the  slope  due  to 
bending  and  the  warping  Y  are  all  constrained  to  vanish,  whereas  the  total  slope  X  and 
the  shear  strains  Y  and  Yf  are  non-zero. 

0  X 

An  element  stiffness  has  also  been  developed  for  a  quadratic  variation  of  shear  strain  by 
adding  two  degrees  of  freedom  in  shear  at  an  internal  node  given  by  £  =  1/2.  The  resulting 
10  x  10  matrix  is  derived  and  transformed  in  a  fashion  exactly  analagous  to  that  for  the  linear 
shear  element.  In  the  quadratic  case  the  static  condensation  is  applied  to  the  four  “internal” 
degrees  of  freedom  so  the  stiffness  employed  in  the  assembly  process  remains  6x6. 
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SECTION  m 

AXISYMMETRIC  SANDWICH  SHELL  ELEMENT 


The  doubly  curved  element  developed  by  Khojasteh-Bakht  (Reference  17)  is  shown  in 
Figure  3.  Displacements  in  the  local  rectilinear  coordinate  system  £  -  77  are  designated  by 
u1  and  u2.  By  matching  the  position,  slope  and  curvature  at  the  nodes,  the  geometry  shown 
in  the  figure  is  substituted  for  that  of  a  segment  of  an  arbitrary  rotational  shell.  Hence,  the 
meridian  of  the  substitute  element  is  defined  by  the  curve 


V  =  £ ( I  “  £>  (a,+  o2  £  +  a3£  Z-I-d4  £3) 


(12) 


where  £  is  a  normalized  coordinate  which  takes  the  values  0  and  1  at  nodes  i  and  j,  re¬ 
spectively.  The  values  of  the  constants  are  given  by 

a  ,  =  tan  . 

a2  =  ton /3.  +-  77"  /2 

a  ,  =  3  (tan  j8  +  t  on  (3  )  —  (77"  -  77"  )  /2 
i  j  J  i 

a3  =  -(5  ton  /3.  +  4  ton  /3.  )  +•  (77."/ 2  —  77") 

_«  .  . _ i _ 

d  £z  R,  cos3/3 

The  displacement  transformation  equations  are  also  necessary: 


u  =  u (  cos  /3  +  u2  sin  ft 

n  .  (13) 

w  =  u(  sin  p  -'ll  cos  p 

A  sandwich  element  is  depicted  in  Figure  4.  Under  axisymmetric  loading,  the  rotation 
of  the  tangent  to  the  meridian  is  given  by 


X  = 


dw _ u 

d  s  R| 


(14) 
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Where  is  the  radius  of  meridional  curvature.  As  in  the  beam  analysis,  this  rotation  can  be 
considered  to  have  two  components.  The  pure  shear  contribution  at  any  location  along  the 
meridian  is  the  same  as  that  given  in  Equation  2  for  the  beam: 

x  s  1  4r" =  x  -  xb  1  rc  hc/d^/f  hf  /d  (15: 

If  u  and  w  are  the  displacements  of  the  reference  surface  and  if  £c  and  C,  ^  are  the 
thickness  coordinates  with  origins  at  the  respective  middle  surfaces  of  the  layers,  the 
kinematic  assumptions  result  in  the  following  expressions  for  the  shell  displacements  due  to 
axisymmetric  loading: 

=  y  ~£c(X-  Yc  ] 

U;-b  =  u  -£f(XTyf)lfx4yc?-2-  (16 

wc  =  w<f  =  wf  =  w 


The  shell  strains  may  be  written  in  the  following  form 


(17) 


where  /3  represents  either  s  or  Q  ,  and  k  either  c  or  f.  Therefore,  applying  the  strain- 
displacement  equations  from  the  classic  linear  thin  shell  theory  (Reference  19)  to  Equa¬ 
tions  16,  one  obtains  the  strain-displacement  equations  for  the  sandwich  shell 


0  jJu _ w_ 

€sc  ds  Rj 

e°  s  ~~  (  u  cos  <f>  -  w  sin  ) 

9c 

dX  ,  dYc 

sc  ds  ds 


COS 

*  '  r  !  x  -  yc  > 

Qt,b  du  w  +  d  dx  -  hc 

dxc 

?  h. 

d/f 

€sf  ds  R,  ”  2  ds  2 

ds 

2 

ds 

t  b 

e°sf’  =  -y- 1  u  cos  <f>  —  w  sin  ) 

+  ■ 

cos  ^ 

d 

IX’ 

Xc  " 


2 


) 
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ds  ds 


(  18) 


However,  in  order  to  employ  the  substitute  element,  the  strains  must  be  expressed  in  terms 
of  uT  and  u2,  the  displacements  in  local  rectilinear  coordinates.  This  can  be  accomplished  by 
substituting  the  transformations  of  Equation  13  into  Equations  14  and  18.  The  lengthy  equations 
that  result  will  not  be  given  here  (Reference  22). 


In  applying  a  stiffness  analysis,  the  assumed  displacement  patterns  are  taken  in  the  local 
rectilinear  coordinate  system  and  are  expressed  in  terms  of  generalized  coordinates.  For  a 
linear  variation  of  shear  strains  and  meridional  displacement  and  a  cubic  variation  of  trans¬ 
verse  displacement,  the  displacement  model  is  given  by 

ut  =  a,  +-a2£ 

Uz  1  a3  +a4^+‘°t5  £2+a6£3 

o  <  £  <  I  (19) 

yc  =  a7  +  a^ 

y, 1  a,  +aio{ 

which  may  also  be  written 

u  (£)  =  <*><£)<* 

5X1  5  X  10  10X1 

uT  =  <  u.  u 2  xb  r  xf  >  (20) 

By  applying  the  strain-displacement  equations  to  Equations  19  and  20,  a  result  analagous 
to  Equation  8  for  the  beam  element  is  obtained 

«  =  B(£)  a  (21) 

(5  X  I  15  X  10  IOX  I 
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In  addition,  it  is  possible  to  construct  the  sparse  matrices  C  and  Z  £  for  the  shell  analagous 
to  the  matrices  in  Equations  6  and  9,  respectively. 

a-  --  c*  (22) 

•  =  Z( C  >« 

The  constitutive  matrix  C  follows  from  the  assumption  of  a  state  of  generalized  plane 
stress.  The  thickness-geometry  matrix  Z  is  derived  directly  from  Equation  17.  Therefore, 
the  element  stiffness  matrix  in  generalized  coordinates  may  be  obtained  from  the  strain 
energy  expression  and  is  as  follows: 

h/g  r  ( £) 

ka  =2^/'b(J)T/  Z(pT  C  Z  (pdC  8(J>  <•?  1231 

10  X  10  O  -h/2 

The  inner  integral  is  readily  evaluated  in  closed  form.  However,  the  integral  with  respect  to 
£  is  more  practically  calculated  by  numerical  means  such  as  Gauss’  integration  formula. 

The  stiffness  matrix  must  be  transformed  to  a  global  coordinate  system  for  the  assembly 
process.  In  addition,  interelement  compatibility  needs  to  be  preserved  only  on  the  displace¬ 
ment  parameters  u,  w,  Xb  and  X.  Hence,  the  element  displacements  in  global  coordinates  are 
selected  as  follows: 

fT  =  <  ui  wi  Xbi  Xi  uj  wj  xbj  Xj  |  Xfi  Xfj  > 


The  transformation  is  given  by 


k  1  T  *a  T 

iox  io 


A'1  T 


Here  A  is  the  transformation  relating  the  local  rectilinear  coordinates  to  the  generalized 


coordinates 


<X>  (  0) 
<X>  (  I  ) 


a  =  a  a 


and  T  is  the  simple  point  transformation  connecting  the  local  and  global  systems 

q  «  T  r 
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Static  condensation  is  employed  to  obtain  the  8x8  stiffness  used  in  the  direct  stiffness 
method.  It  should  be  noted  that  if  the  shell  has  discontinuities  in  the  meridional  slope,  the 
curvilinear  displacements  u  and  w  are  not  uniquely  defined  at  these  points.  In  this  case,  the 
translations  in  global  coordinates  can  be  taken  in  the  rectilinear  directions,  i.e.,  u  and  u 
in  Figure  4.  r 

For  rotational  shells  which  are  closed  at  the  end,  a  special  cap  element  is  used  to 
represent  the  closure.  The  stiffness  for  this  element  is  obtained  in  basically  the  same  manner 
as  for  the  frustum  element.  However,  “internal  boundary  conditions”  are  applied  at  the  apex 
to  reduce  effectively  the  number  of  degrees  of  freedom  (Reference  17).  From  the  requirement 
of  axisymmetry  it  is  apparent  that  these  conditions  at  the  apex  node  are : 

ur  :x  :/c:n:°  01  r 

In  addition  to  the  linear  shear  model,  element  stiffnesses  have  also  been  derived  for  a 
quadratic  variation  of  shear.  As  for  the  beam,  an  internal  node  at  £  =  1/2  is  utilized  and 
static  condensation  enables  reduction  of  the  stiffness  from  12  x  12  to  8  x  8. 
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SECTION  IV 

LUMPED  MASSES  FOR  DYNAMIC  ANALYSIS 


The  mass  of  the  structure  is  concentrated  at  the  nodal  points  so  that  the  inertial 
properties  can  be  represented  by  a  diagonal  mass  matrix.  Felippa  (References  20  and  11) 
has  demonstrated  that  this  lumped  mass  procedure  provides  satisfactory  fundamental 
frequencies  and  mode  shapes  with  less  computational  effort  than  the  consistent  (or  distributed) 
mass  approch.  For  the  same  number  of  nodal  points,  the  former  method  results  in  fewer 
equations  for  the  eigenvalue  problem.  When  meshes  are  arranged  so  there  are  the  same 
number  of  eigenproblem  equations  for  both  techniques,  Felippa’ s  results  indicate  that  the 
lumped  mass  approach  produces  the  more  accurate  frequencies. 

Rotatory  inertia  has  been  shown  to  be  a  more  important  factor  in  the  vibration  of  sandwich 
plates  (Reference  8)  than  in  the  dynamics  of  homogeneous  plates.  For  translational  displace¬ 
ments  it  is  sufficient  to  idealize  the  lumped  mass  as  a  point  on  the  reference  surface 
(Figure  5a).  However,  for  rotatory  effects,  the  distribution  of  the  mass  through  the  depth  of 
the  beam  must  be  maintained.  This  is  especially  true  for  sandwich  structures  where  the 
outermost  layers,  the  facings,  may  be  much  denser  than  the  core.  Hence,  one  can  visualize 
the  mass  as  being  lumped  along  the  material  line  originally  normal  to  the  reference  surface, 
with  no  concentration  of  the  mass  across  the  depth  (Figure  5b).  In  effect,  this  is  the  same  as 
multiplying  the  mass  moment  of  inertia  of  the  cross  section  by  the  tributary  area.  It  should 
be  noted  that  the  rotatory  inertia  is  associated  with  the  rotation  of  the  normal  to  the  reference 
surface,  i.e.,  X^. 
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NODE  i-l  NODE  i  NODE  i  +  l 

j  i  i 

_ i _ wmfflMm/m  1 — 
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1  '  ’ 

1  i 

_  TRIBUTARY  li/2  I 

1  LENGTH  i  i 

(a)  LUMPED  MASS  FOR  TRANSLATIONS 


(b)  LUMPED  MASS  FOR  ROTATION 


Figure  5.  Arrangement  of  Lumped  Masses 
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SECTION  V 
EXAMPLES 

The  above  finite  element  formulation  has  been  applied  to  various  sandwich  beam  and 
shell  problems  and  the  results  compared  to  solutions  from  other  methods  and  sandwich 
theories.*  A  sampling  of  these  problems  is  presented  in  this  section  to  demonstrate  the 
efficacy  of  the  method.  In  general,  the  displacements,  stress  resultants  and  natural  fre¬ 
quencies  from  the  finite  element  method  compare  favorably  to  corresponding  quantities 
obtained  by  established  theories. 

END-LOADED  CANTILEVER  BEAM 

A  cantilever  sandwich  beam  of  unit  width  with  a  unit  load  at  the  free  end  illustrates  the 
effect  of  a  constraint  on  the  warping.  The  dimensions  and  properties  are  selected  as  follows: 

hc  =  0.  5"  ,  h  f  =  0  04"  ,  h  =  0.  58" 

Ef  =  I07  psi  ,  Gf  =  4  X  10®  psi ,  Xf  =1 
Ec  =  2  X  I04  psi  ,  Gc  =  I04  psi,  \c  =  I 
span  L  *  10",  load  P  =  1.0  lb. 

Evenly  spaced  meshes  of  5  and  10  elements  as  well  as  uneven  meshes  are  used  for  both  linear 
shear  strain  (L  elements)  and  quadratic  shear  strain  (Q  elements). 

The  displacement  solution  for  5-L  elements  is  shown  in  Figure  6;  the  displacements 
from  a  5-Q  analysis  fall  about  midway  between  the  5-L  result  and  the  solution  of  Reference  7. 
For  all  meshes  and  elements  used,  the  overall  stress  resultants  are  correct  to  about  five 
significant  figures,  so  these  results  are  not  shown  graphically.  Of  interest,  however,  is  the 
distribution  of  shear  force  between  the  facings  and  core.  The  fraction  of  the  shear  assumed 
by  the  core  is  shown  in  Figure  7.  The  theory  from  Section  1.2  of  Reference  6  does  not  take 
into  account  either  the  warping  behavior  or  the  bending  stiffness  of  the  facings;  it  assumes 
that  all  shear  is  taken  by  the  core.  The  refinement  to  take  into  account  the  restraint  on 
warping  and  the  consequent  flexure  of  the  facings  about  their  own  middle  surfaces  is  given 


*For  distributed  loadings,  consistent  nodal  loads  are  employed;  e.  g.,  for  the  beam, 
q  -  j  1  { £  )T  d£  ,  where  Q  are  the  nodal  loads  corresponding  to  displacements 

q  and  “p  are  the  specified  loads. 
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Figure  6.  Deflection  of  a  Cantilever  Beam 


Figure  7.  Shear  Distribution  Between  Face  and  Core 
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in  Section  1.3  of  Reference  6.  This  formulation  is  due  to  van  der  Neut.  Finally,  Yu’s  theory 
(Reference  7)  considers  both  the  warping  and  the  shearing  of  the  facings  and  thus  gives  the 
distribution  of  shear  among  the  various  layers.  Figure  7  demonstrates  that  with  a  proper  mesh 
refinement,  the  finite  element  method  gives  an  adequate  representation  of  this  phenomenon. 
Moreover,  the  quadratic  shear-strain  elements  enable  a  satisfactory  representation  with 
fewer  elements.  At  the  support,  the  constraint  against  warping  causes  most  of  the  shear  to  be 
carried  by  the  facings. 

Failures  have  been  found  to  occur  in  the  facings  near  fixed  supports  of  aerospace 
sandwich  structures.  For  this  reason,  the  facing  layers  are  usually  doubled  in  thickness  in 
these  regions.  The  above  results  give  an  insight  into  the  shear  redistribution  which  neces¬ 
sitates  the  use  of  such  doubler  plates.  In  fact,  the  finite  element  method  is  suited  for  design 
of  doubler  plates  since  the  computer  program  is  readily  modified  to  account  for  elements  with 
differing  face  thicknesses.  Hence  it  is  possible  to  include  these  reinforcing  layers  in  the 
analysis. 

HEMISPHERICAL  SHELL  UNDER  MEMBRANE  LOAD 

To  check  the  effectiveness  of  the  basic  element  and  of  the  computer  program  for  shells, 
the  membrane  states  of  both  cylindrical  and  spherical  shells  have  been  investigated.  Generally, 
the  results  are  satisfactory  in  that  both  deflections  and  stress  resultants  agree  with  theoret¬ 
ical  values.  A  typical  example  is  presented  here.  The  sandwich  hemisphere  has  the  following 
properties: 

hc  =  0 . 5",  h  f  =0,  04" 

E f  =  I07  Psi,  vf  *  0.3,  Gf  =  3.85  X  10®  psi,\f  =  1.0 
Ec  =  2.  6  X  !04  psi  ,  vz  -  0.  3  ,  Gc  =  I04  psi  ,  Xc  =  1.0 
radius  a  =  100",  load  pz  =  -1.0  psi 

Three-  and  nine-element  representations  are  used  withbothlinear  shear  and  quadratic  shear. 
Results  are  essentially  the  same  for  the  two  shear-strain  models,  so  only  the  solution  using 
the  less  refined  model  is  presented  here. 

When  roller  supports  that  restrict  only  the  meridional  displacements  at  the  free  edges  are 
used,  it  is  easily  shown  that  the  theoretical  solution  is 

N  s  =  Ng  =50  lb  /  in. 
w  =  -  0.004305  in. 
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The  finite  element  solution  for  three  elements  is  shown  in  Figure  8.  It  is  seen  that  the  results 
agree  very  closely  with  the  theoretical  answers.  The  nine-element  solution  is  even  better, 
and  is  not  shown  since  it  does  not  differ  significantly  from  the  exact  solution. 


SHALLOW  SPHERICAL  CAP  WITH  PARTIAL  DISTRIBUTED  LOADING 

The  last  static  shell  problem  presented  here  is  a  simply  supported  shallow  spherical  cap 
subject  to  distributed  loading  over  a  portion  of  its  surface.  The  closed-form  solution  for 
this  problem  in  terms  of  Thomson  functions  has  been  given  by  Rossettos  (Reference  21).  He 
neglects  the  shearing  of  the  faces  and  the  bending  and  extension  of  the  core,  so  the  following 
properties  are  selected: 


=  0.  95"  , 

h  f  =  0.02 

5"  ,  h  =  1  " 

.si  ,  vf  = 

03,  Gf 

=  I020  psi,  \f  =  1 

:  =  0  ■  GC 

=  I05  psi  , 

V 

□  -  20  , 

su  ppor  t  ed 

edge  at  c p  =15° 

A  uniform  load  of  1  psi  is  applied  in  the  axial  direction  over  that  portion  of  the  surface  given 
by  0<  <£<3  degrees. 


The  cap  is  analyzed  using  5  and  10  linear-shear  elements.  Deflection  results  are  presented 
in  Figures  9  and  10,  and  stress  resultants  are  plotted  in  Figures  11  through  13.  In  general, 
satisfactory  results  are  obtained  from  the  five-element  representation.  The  small  difference 
between  the  five-  and  ten-element  results  indicates  that  the  finite  element  solution  has 
effectively  converged. 

No  examples  are  presented  for  shells  which  are  not  shallow  because  of  the  absence  of 
published  solutions  with  which  to  compare  them,  although  the  computer  program  developed 
is  capable  of  solving  such  problems. 


SIMPLY  SUPPORTED  CYLINDRICAL  SHELL 

One  type  of  shell  for  which  solutions  for  natural  frequencies  are  readily  available  is  the 
cylinder.  Yu  (Reference  9)  has  derived  a  three-branch  frequency  equation  for  an  infinite 
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Figure  8.  Hemispherical  Shell  Under  Membrane  Load  (Three  Elements) 


ANGLE  <j>  (  DEGREES) 


Figure  9.  Radial  Dispacement  of  a  Shallow  Sphere 


ANGLE  <p  (  DEGREES  ) 
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Figure  10.  Meridional  Displacement  of  a  Shallow  Sphere 


ANGLE  <j>  (  DEGREES  ) 


Figure  11.  Tangential  Stress  Resultants  of  a  Shallow  Sphere 


ANGLE  <f)  (DEGREES) 
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Figure  12.  Shear  Stress  Resultant  of  a  Shallow  Sphere 
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Figure  13.  Moments  of  a  Shallow  Sphere 
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cylinder  in  an  extension  of  his  theory  for  sandwich  plates  (References  7  and  8).  This  equation 
is  also  applicable  to  a  simply  supported  shell  of  finite  length.  A  cylinder  with  the  following 
properties  is  now  analyzed. 

h  =  0.5"  ,  h  --  0.025“  ,  h  =  0.55" 
c  ’  f  ' 

Ef  =  I07  P  5  i  ,  Vf  -  0.  3  ,  Gf  =  3.8  5  X  I06  psi,  X  =  I 

Ec  =  2.6  X  )04  psi,  =  0.  3  ,  Gc  =  I04  psi,  Xc  =1 

p  -  0.1  lb/  in.3  ,  p  =0.005  lb  /  in.3 
rf  rc 

radius,  a=  20"  ,  span,  L  =  10 

As  in  the  preceding  examples,  these  properties  are  typical  of  sandwich  construction. 

The  natural  frequencies  computed  by  the  finite  element  method  are  compared  with 
solutions  by  Yu  in  Table  I.  Rotatory  inertia  is  included  in  all  cases  and  quadratic -shear 
elements  are  used.  Many  more  frequencies  than  would  be  of  practical  interest  are  shown  in 
the  table  in  order  to  evaluate  better  the  overall  effectiveness  of  the  finite  element  approach. 
It  is  noteworthy  that  the  lowest  frequencies  of  each  of  the  three  modes  are  approximated 
regardless  of  the  number  of  elements  used.  The  number  of  frequencies  given  by  the  finite 
element  method  for  each  mode  depends  on  the  number  of  degrees  of  freedom  available  of  the 
type  that  are  necessary  to  characterize  the  particular  mode. 

For  shell  structures  there  is  usually  more  than  one  branch  of  the  frequency  equation  that 
is  of  engineering  importance.  Careful  study  of  the  finite  element  mode  shapes  must  be  under¬ 
taken  in  order  to  identify  the  frequencies  with  the  appropriate  branch.  This  is  especially  true 
in  preliminary  analyses  wherein  the  frequencies  have  not  yet  converged  to  a  predictable 
pattern. 

In  the  present  example,  simple  supports  which  preclude  translation  in  any  direction  are 
used.  Hence  the  breathing  mode  or  fundamental  radial  expansion  mode  is  prevented.  However, 
the  example  has  been  recomputed  with  supports  that  restrain  only  longitudinal  displacements 
in  order  to  obtain  an  estimate  of  the  cut-off  frequency  of  this  radial  mode.  Yu’s  solution 
(Reference  9)  for  this  frequency  is  8520  radians/second.  The  finite  element  approximations 
for  five,  ten,  and  twenty  elements  are  8360,  8500  and  8520  radians/second,  respectively. 
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TABLE  I  —NATURAL  FREQUENCIES  OF  A  SIMPLY  SUPPORTED 


SANDWICH  CYLINDER  (RAD. /SEC.) 


MODE 

TYPE 

R 

Finite  Element  Method 

5  elems . 

10  elems . 

20  elems . 

1 

LI 

8740 

9040 

9080 

9080 

L2 

11900 

11900 

12100 

12200 

L3 

16300 

15500 

16700 

16800 

l4 

21000 

18000 

21500 

21900 

L5 

25900 

26000 

27000 

L6 

30800 

30000 

32300 

L7 

35700 

33400 

37500 

L8 

40700 

35900 

42700 

L9 

45700 

37400 

47700 

L10 

50700 

52700 

R1 

53600 

52600 

53300 

53500 

Lll 

55700 

57400 

13 

L12 

60700 

61900 

Ik 

L13 

658OO 

66000 

15 

TSCO 

69800 

71500 

72800 

73100 

16 

LI  k 

70800 

698OO 

IT 

L15 

75800 

73100 

18 

LI  6 

80800 

75900 

19 

LIT 

85900 

78100 

20 

LI  8 

90900 

79700 

21 

TS1 

92500 

91200 

94100 

94900 

22 

L19 

95900 

80700 

23 

L20 

101000 

2k 

L21 

106000 

25 

R2 

107000 

99800 

105000 

107000 

32 

TS2 

140000 

128000 

138000 

141000 

37 

R3 

161000 

137000 

155000 

159000 

1+7 

TS3 

195000 

187000 

194000 

R4 

214000 

200000 

211000 

TS4 

252000 

235000 

249000 

R5 

268000 

241000 

261000 

TS5 

311000 

278000 

303000 

r6 

321000 

275000 

309000 

TST 

370000 

356000 

R7 

375000 

_ 

356000 

L  =  Longitudinal,  R  =  Radial,  TS  =  Thickness  Shear 
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SECTION  VI 
CONCLUSIONS 


The  finite  element  method  has  been  extended  to  the  refined  analysis  of  multilayer  beams 
and  axisymmetric  shells.  In  the  theory  employed  no  restriction  is  placed  upon  the  ratios  of 
the  layer  thicknesses  and  properties.  The  method  is  applicable  to  structures  wherein  shearing 
deformations  are  significant,  including  sandwich  construction. 

One  of  the  features  of  the  formulation  developed  herein  is  the  capability  of  approximating 
the  warping  phenomenon.  Hence,  the  distribution  of  the  shearing  force  among  the  various 
layers  can  be  determined.  Another  feature  is  the  use  of  lumped  rotatory  inertia  in  the  dynamic 
analyses.  This  type  of  inertia  is  a  prerequisite  for  the  inclusion  of  the  thickness-shear  modes 
of  behavior,  which  are  important  for  some  soft-core  sandwich  structures.  A  third  aspect  of 
the  present  work  is  the  possibility  of  neglecting  either  the  shearing,  extension  or  bending 
of  any  individual  layer.  By  use  of  this  capability,  the  effects  of  various  approximations  can  be 
evaluated  for  different  geometrical  or  material  properties. 

Generally,  other  available  analysis  techniques  can  only  be  applied  to  sandwich  structures 
with  the  simplest  configurations.  Hence  there  is  a  great  potential  for  the  application  of  the 
finite  element  method  to  the  solution  of  sandwiches  of  arbitrary  shape. 
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